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ABSTRACT 

We present a procedure for determination of positions and orbital elements, and 
associated uncertainties, of outer Solar System planets. The orbit-fitting procedure 
is greatly streamlined compared to traditional methods because acceleration can be 
treated as a perturbation to the inertial motion of the body. These techniques are 
immediately applicable to Kuiper Belt Objects, for which recovery observations are 
costly. Our methods produce positional estimates and uncertainty ellipses even in the 
face of the substantial degeneracies of short-arc orbit fits; the sole a priori assumption 
is that the orbit should be bound or nearly so. We use these orbit-fitting techniques 
to derive a strategy for determining Kuiper Belt orbits with a minimal number of 
observations. 

Subject headings: Kuiper Belt — celestial mechanics 

1. Introduction 

The discovery of the Kuiper Belt (Jewitt and Luu 1993) created a new vista upon the formation 
and early evolution of the Solar System. Study of these objects' basic physical properties of course 
requires sufficiently well-determined orbital parameters to deduce a future position and the distance 
to the object. The orbital parameters are themselves of great interest, with the phase space of orbits 
threaded by an intricate web of resonant and long-lived chaotic zones (Malhotra 1996). For both 
the practical goal of object retrieval and the higher goal of understanding the dynamical structure 
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of the population, it is essential to not only have accurate orbits but also to quantify the uncertainty 
in the orbit. 

There are centuries-old methods for determination of orbital parameters from observed posi- 
tions. Present-day desktop computers execute these solutions trivially, and a linearized propagation 
of the positional uncertainties into the orbital elements is straightforward as well (Muinonen and 
Bowell 1993). Workstations are even fast enough to bound the uncertainty region with a brute-force 
sampling of the 6-dimensional orbit space in many cases. Why then should we bother developing 
another orbit-fitting technique? KBOs pose a particular challenge because the objects are faint, so 
recovery observations are quite costly, requiring substantial investment of 2-4-meter telescope time 
for all but the brightest known objects. As a consequence, most known objects have been observed 
only a few times, leading to substantial degeneracies in the orbit fits. We need to know the error 
ellipse (in sky position and in orbital parameter space) even in the face of these degeneracies. In 
the best case a brute-force method will work on underdetermined orbits but leave us ignorant of 
the nature of the degeneracy. In the worst case the degeneracies will lead to long or effectively 
infinite computation times. 

The Minor Planet Center (MPC) (http://cfa-www.harvard.edu/iau/mpc.html) provides or- 
bital elements and predicted ephemerides for new objects. The approach of the MPC to short-arc 

objects is to select the simplest "sensible" orbit that fits the data — e.g. a circular orbit, if it fits 
the data and docs not imply a close approach to Neptune; a 3:2 resonant orbit with Neptune if the 
circular orbit would imply a close approach; or, if the circular orbit is a poor fit, a "Vaisala" orbit, 
with the object near perihelion (Marsden 1985, 1991), if this does not imply a close encounter. The 
MPC orbits in effect remove the degeneracies by assuming that the orbit is most likely to resemble 
those of known objects, and consequently predict positions that in most cases are quite close to the 
actual recovery positions. There are a few difficulties with this approach, however: first, focussing 
recovery efforts on the positions predicted by these favored orbits will inevitably bias the recovered 
population toward such orbits. This makes it difficult to quantify the portion of the population 
which may be in unusual orbits as they are preferentially lost. Second, uncertainties on positions 
and orbital elements are not currently available from the Minor Planet Center. Positional uncer- 
tainties are essential for proper planning and cvahiation of recovery observations. Uncertainties on 
orbital elements are important for analyses of the dynamical characteristics of the population. For 
example, there was speculation that the 2:1 resonance with Neptune was empty or underpopulated, 
but 1997 SZio and 1996 TRgg were reclassified from 3:2 resonators to 2:1 resonators after further 
observations. Finally, it is sometimes necessary to calculate orbits "on the fly" during an observing 
run, and the MPC should not be expected to provide instantaneous analyses. 
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2. Methods 
2.1. Motivation 

Our approach to orbit-fitting for KBOs is motivated and guided by the foUowing differences 
from the historically more common application to the asteroid population: 

• Recovery observations are very expensive, as noted. This leads us to investigate, in §4, the 
minimal number of observations required to reach a given accuracy in orbital elements. The 
prediction of future positions should be a stable procedure even in the presence of nearly- 
degenerate orbital elements. This will require that we fully understand the nature of the 
degeneracies from short arcs. Wc would like to find a parameterization of the orbits in which 
the degeneracies are confined to as few parameters as possible. 

• The release of the USNO-A2.0 astrometric catalog^ now makes it possible to measure astro- 
metric positions to (y.'2 accuracy over the full sky (Dcutsch 1999) in a typical CCD image, and 
to produce reasonable estimates of the uncertainties on each position measurement. Uncer- 
tainty estimation can therefore proceed by the straightforward propagation of measurement 
errors. Space telescopes and adaptive-optics ground-based telescopes will commonly produce 
relative positions to Of.'Ol accuracy in the near future. Very accurate positions should in 
principle produce accurate orbits even over relatively short arcs. 

• KBOs are at distances d > 30 AU so their apparent motions are dominated by reflex motion. 
The observed arcs are a small fraction of the orbital period P even after a decade. Apparent 
motion patterns are thus very simple and unambiguous in comparison, say, to Near-Earth 
Objects. 

• If we parameterize the distance to the object by 7 = (1 A\J)/d < 0.03, the acceleration of the 
KBO is 7^ < 10~^ times smaller than Earth's acceleration, and the component transverse 
to the line of sight is 7^ ~ 10^^'^ times Earth's. The KBO motion is nearly incrtial, and 
gravitational acceleration can be treated as a perturbation. Instead of parameterizing orbits 
by the usual element vector a = {a,e,i,Q,u>,Tp}, an orbit is more stably specified by some 
Cartesian phase space vector P = {xq, yo, zo,xo,yo, zq} at the time of initial observation. 

• The annual parallax is limited to 27 « 3°, and the total object motion is only ^ 10° over 
a decade. We will therefore express sky positions as coordinates {Ox,Oy) in a tangent-plane 
projection of the celestial sphere about an appropriate reference point, e.g. the first observed 
location. Furthermore, the line of sight to the object changes by only a few degrees, and we 
will find it convenient to align our Cartesian coordinates with the z axis along the initial line 
of sight. 



^The USNO-A2.0 catalog, produced by D. Monet et al, is described at 
http://www.nofs.navy.mil/projects/pmm/USNOSA2doc.html. 
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Fig. 1. — Panel (a) illustrates the relation of our angular coordinate system to ecliptic coordinates. Our 
coordinates are a tangent-plane projection of the sky about the initial line of sight, (^Oi^o)- Panel (b) 
illustrates the spatial coordinate system we adopt in relation to the ecliptic spatial system Xoc. Our 
system has origin at the observer location on Earth's surface for t = and the z axis along the initial line of 
sight to the target. The vectors and xb toward the observer and solar system barycenter are illustrated 
as well. Both spatial and angular x and y axes are aligned to local ecliptic East and North, respectively. 

• Daily parallax is < 8 x 10^^7 < (X'S. For most ground-based observations this will be too 
small to yield useful information. 

2.2. Exact Equations 

The input data are the observed right ascension and declination (ai,5i) of the target at times 
ti. We will assume that all positions have been measured in the J2000 celestial reference frame 
via astrometry tied to background stars in the USN0-A2 catalog. Measuring relative to the back- 
ground stars exactly cancels the effects of stellar aberration, and reduces the maximal error from 
gravitational deflection of starlight to the milli-arcsecond level, so we can treat the problem with 
Euclidean geometry. We transform the observed coordinates to a tangent-plane projection {Ox,Oy) 
about a reference location (ao,5o)- The 6x axis is taken in the Eastward direction of the J2000 
ecliptic coordinate system, and By points toward J2000 ecliptic North. Formulae for the conversion 
from (a, (5) to {Ox,Oy) are straightforward and summarized in Appendix A. Figure 1 illustrates the 
relation of our coordinate systems to the ecliptic. 

We will generally take this projection axis to be the first observed position, and will take the 
first observation to be at time to = 0. We will assume that 9x and 9y have equal and uncorrelated 
measurement uncertainties cjj (generalization to non-circular error regions is straightforward). 

We define our inertial Cartesian coordinate system by placing the origin at (or near) the 
location of the observer at the initial epoch t = 0. The z axis is directed toward the reference 
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direction (ao,5o) of the tangent-plane projection, i.e. along the initial line of sight. The y axis is 

in the plane containing the z axis and the J2000 North ecliptic pole, i.e. parallel to the 6y axis of 
the tangent plane. The x axis then lies in the ecliptic plane and points along 9^ to ecliptic East. 
Unless otherwise specified, we will take the units of distance and time to be AU and Julian year, 
respectively. The coordinate systems are illustrated in Figure 1. 

If x(i) represents the position of the target at time t and X£;(t) is the position of the observer 
on Earth, then the observed sky position at time t is 



^ ^^-j _ x{t-At)-XE{t) 

^y(^) = z(t-At)-ZE{t) 



z{t-At)-ZE{t) (-,\ 
y{t-At)-yE{t) 



where At is the light-travel time |x(t — At) — XB(t)|/c. The definition of At is recursive but in 
practice is easily solved for iteratively because of the lack of ambiguity in KBO orbit solutions. 

We wish our uncertainties to be limited by the errors in the astrometric catalog, so our re- 
quirement is to model the target position to < (X'l accuracy over k. 10-year time spans. For target 
distance > 30 AU, the object and observer positions must be correct to ~ 2000 km. The Earth's 
position relative to the Sun and Solar System barycenter (SSB) is available from the JPL DE405 
Ephemeris^ (Standish et al. 1995) to an accuracy many orders of magnitude better than we require. 
A simplistic model of the Earth as an oblate spheroid in constant rotation gives the observatory 
position to accuracy well within our tolerance, so we may consider xg to be known exactly. The 
interpolated JPL ephemeris and topocentric correction are rapidly calculable. 



2.2.1. Orbit Approximations 

We write the target orbit as 

x(t) =xo + xot + g(t), (2) 
where g is the gravitational perturbation satisfying 

g(t = 0) = g(t = 0) = 0, (3) 
_ ^ -GM,[x(t)-x,(t)] 
- ^ |x(t)-x,(t)|3 ' 

where the sum runs over the other bodies in the Solar System. We include forces from the Sun 
and the giant planets, obtaining their positions from the DE405 ephemeris. For objects beyond 
30 AU, one could approximate the total gravitational field of the Solar System as a monopole 
originating the Solar System barycenter, with errors of only ~ 10~^ AU over 10 years (barring 
Neptune encounters), yielding positions accurate to < 1". The computational expense of the full 



^The DE405 ephemeris is available at http://ssd.jpl.nasa.gov/eph_info.html 
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N-body form Eq. (4) is insignificant, however, as the time steps in the orbit integration can be 
quite long (20 days or more) without significant loss of accuracy over a decade. 

The combination of Eqs. (1), (2), and (4) define our model for the position {0x,9y) as a 
function of time and the initial state vector P. A comparison of our positions for Pluto, Neptune, 
and 1992 QBl to those generated by the JPL ephemeris indicate residual differences of < C/.'OS for 
a 10-year integration. 



2.3. Approximate Formulae &; Degeneracies 

We can understand the nature of the short-arc degeneracies by expanding Eq. (1) in powers 
of time t and distance parameter 7. The formulae are simpler if we redefine our phase space basis 
in terms of the vector a = {a, (3, 7, a, $, 7}, where 

a = xo/zo f3 = yo/zo 7 = 1/^0 .^x 
a = xo/zo $ = yo/zo ^ = zo/zo. 

Please note that a, /5, and 7 quantify the motion along our three axes but are not the time 
derivatives of a, (3, and 7. We have chosen our coordinate system so that a, /? ~ 0, and we expect 
d,/3,7 ~ 27r7^/^ for bound orbits. If the orbit is nearly circular, the line-of-sight velocity is down 
another factor of 7 so that 7 ~ I-k^^I"^ . The x component of the apparent motion is now 

l + 7t + 75^(t)-7ZE(i) ^' 
K, a -I- (a - ^XE)t (7) 

+7(5x - Xe)^ I'i' - (7 - lZE)[pL - JXE)t^ 

In the second line we have made the approximation that t <^ 1 year. The equation for 9y is 
analogous. We see that the solution for the orbital parameters a has three regimes: 



2.3.1. Slope Regime 

For arcs spanning consecutive nights, we can only determine the instantaneous positions a and 
P, plus the slopes 

9x = {a- JXe) ~ a - 7(27r cos ^0) (8) 
Oy = 0-7yE)^P (9) 

The approximations assume that we are observing near the ecliptic plane, and that the Earth's 
orbit is circular with phase ^0 at time t = relative to the Sun-target vector. In this situation 
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the line-of-sight motion 7 is completely unconstrained, and there is a total degeneracy between the 
distance 7 and the ecliptic motion a. For a circular orbit, \6l\ ~ (^/7/ cos (/)o)(7i;£'), and for any 
bound orbit |d| < (-^27/ cos 0o)(7^£;)- The usual observational strategy, therefore, is to observe 
near opposition (^0 = 0) to minimize the proper-motion term a. relative to the parallax term 7Xe, 
giving an estimate of 7 to fractional accuracy of 1/27 « 0.25 or better. An ephemeris may be 
predicted under the assumption of a prograde circular orbit. Below we will describe a simple way 
to place an error ellipse on the position without assuming a circular orbit. 



2.3.2. Acceleration Regime 

As the arc length grows, the next bit of information gleaned is the apparent acceleration vector: 

= 7(5a; - a^E) - 2(7 - 7ii/)(a - 7i;£;) (10) 
~ 7(ffa; - 47r^sin0o) + 27 (27r7cos0o)- (11) 

The second line again assumes an ecliptic observation and a circular Earth orbit, and keeps leading- 
order terms only. We see first that the effect of 7 on the acceleration is a factor 7/27r < 0.005 smaller 
than the ^xe term. In this regime, therefore, 7 is still indeterminate. The target's transverse 
acceleration is smaller by a factor 7^ than the Earth's acceleration and can also be ignored at 
this point. The apparent acceleration is then 

Qx ~ -7(47r2sin,^o) (12) 

Wc sec that the apparent acceleration is a robust and simple way to constrain 7, the distance to 
the target. This method fails however, for observations at opposition (^0 = 0), where the reflex 
acceleration vanishes. 

What time span is needed to measure the acceleration to a useful accuracy? Let us presume 
that the object is initially detected in a pair of observations separated by 24 hours, then recovered 
in a single observation some time T later, and that each observation has astrometric uncertainty 
a. Using Eq. (12) to determine 7, the fractional error in the distance to the target is then roughly 

^ = (13) 

d 7 V27r27sin(/)o(24h)r 

I ^ \ Aweck\ /sin45°\ ( d \ 

The distance to most KBOs can be determined from the apparent acceleration with only 1 week's 
arc, as long as the observations are not near opposition. For a main belt asteroid, the apparent 
acceleration is an order of magnitude larger and is easily detected in a 24-hour span. This means it 
is easy to distinguish main belt asteroids from KBOs away from opposition, even near the main-belt 
turnaround points. 
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A typical (non-opposition) 1-week arc thus gives a useful estimate of 5 of the 6 phase-space 
parameters of the orbit, with the line-of-sight motion being completely indeterminate. In this basis 
the degeneracy is confined to a single parameter, whereas the degeneracy is shared between most 
of the traditional orbital elements. In particular neither a nor e is at all well constrained at this 
point. 

It is interesting to note that the orbital angular momentum is sensibly constrained in this 
regime, since the line-of-sight velocity is nearly radial to the Sun and contributes little to the 
angular momentum. The energy of the orbit, however, is 

E = -GMqj{1 - 27 cos po + l^Y^''^ + {p? + /?2 + 72)7272, (14) 

and is uncertain until 7 can be constrained. Here /3o is the solar elongation of the target at the 
initial epoch. 



2.3.3. Fully Constrained Regime 

The line-of-sight velocity can be constrained as the Earth moves around its orbit to give us a 
new point of view. From Eq. (6), the leading 7 term in the sky position is 'ytjX}^{t) < 77^. We 
know that 7^ < 2GMq^'^ from the simple constraint that the object be bound. A useful constraint 
on 7 is obtained only after sufficient time that a 7 value below the binding limit would produce 
a measured displacement that is distinguishable from terms in the other 5 parameters. This can 
occur on a timescale of months for high-quality observations. The uncertainty in 7 (or z) dominates 
the phase space uncertainties for essentially all lengths of arc, though it does not always dominate 
the ephemeris uncertainties. 



2.4. Fitting Procedures 

With 3 or more observations spanning several months or more, we will be in the fully- 
constrained regime. At this point the constraint of the orbital parameter vector oc from the mea- 
sured positions 9i becomes a straightforward non-linear minimization problem. The formalism 
for Bayesian orbit fitting is discussed in detail by Muinonen and Bowell (1993). We will assume 
Gaussian measurement errors, in which case the most likely parameter vector a. is the least-squares 
value, minimizing 



X 



The model 9x{ct,ti) for the predicted position is embodied by Eqs. (5), (1), (2), and (4). The 
minimization of can use standard algorithms. We have implemented a slight modification of 
the Levenberg-Marquardt method from Press et al. (1988). For efficient location of the minimum. 
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we require a good starting point for the non-linear search, and estimates of the derivative matrix 
de/d(x. 

The determination of the derivatives is again greatly simplified for distant objects. The gravi- 
tational acceleration is, over a < 10-year period, a small perturbation to the inertial motion, barring 
close encounters. The inertial motion will be ~ 2T\-^t and the acceleration term will integrate to 
~ 47r^7^t^/2. Furthermore the acceleration is toward the barycenter, so the transverse component 
is down by another factor of 7. The angular displacement due to gravity is only a few arcseconds 
per year, > 300 times smaller than the inertial motion and reflex motion. It is clear that in calcu- 
lating the derivative matrix dx/da from Eq. (2) that we can ignore the gravitational term. All of 
the derivatives are thus quite trivial and there is no need to integrate the derivatives of the orbit. 

Selection of a sensible starting point for the minimization is also simple for distant objects. We 
produce a version of the equations of motion which ignores gravity, is linearized, and also ignores 
the 7 term: 

9j; a + at — jXE{t), (16) 

These equations are linear in the parameter vector oc and hence have a closed-form minimization 

requiring only the inversion of a 5 x 5 matrix (Press et al. 1988). These very simple equations are 
actually good to a few arcseconds for the first year of observation, and in any case provide (with 
7 = 0) a good starting point for the non-linear minimization. 

Using the analytic derivatives, the Levcnberg-Marquardt algorithm also calculates the covari- 
ance matrix 5]q of the fitted parameters. Since the model 9{a,t) is dominated by terms linear in 
the parameters, the linearized approximation 

Ax" = (a - &fi:a{oc - d) (17) 

will be accurate even for poorly constrained orbits. This would not have been the case had we 
chosen the orbital-element vector a as the parameter set. 



2.5. Position and Orbital Element Estimation 

Once the best-fit o; vector and its covariance matrix are determined from the fit, determination 
of orbital elements and future positions is also simple. The predicted position at any time t is 9{a, t). 
The uncertainty ellipse for the position is specified by the covariance matrix 

The derivatives, as above, are rapidly calculated by ignoring gravity. Again the position is nearly 
linear in our chosen orbital basis, so this linearized transformation will be accurate even in the case 
of large uncertainties in the parameters. 
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The best-fit orbital element vector a is determined from the best-fit phase-space parameter set 
a via the usual equations, which are reviewed in Appendix B. The determination of an uncertainty 
ellipsoid in orbital-element space requires mapping of the covariance matrix with the derivatives 
da/ da. Some of these derivatives are quite complicated analytic expressions, but can still be coded 
and rapidly evaluated. The map from a to a is non-linear, and in some cases degenerate, so the 
linearized covariance matrix is formally correct only when the uncertainties in the elements become 
small. We do, however, obtain a useful estimate of the element uncertainties even for less exact 
determinations . 



2.5.1. Fitting Singly Degenerate Orbits 

If the arc is too short, the minimization of "will proceed but the best-fit parameters may 
be non-physical. The degeneracy in the "acceleration regime" will be manifested as one or more 
large elements of the covariance matrix. In particular, several-week orbits will have large values of 
S-y-y. We recognize these orbits as unrealistic because they may be unbound. The assumption of a 
bound orbit limits 7 to (c/. Eq. [14]) 

7' < 7bind = GMqj^{1 + 7^-27 cos /Jo)"^/" - - p\ (19) 

A conservative approach to assigning uncertainties in this case is to assume that 7 has a uniform 
a priori probability in the interval [— 7bind) +7bind]- If) however, we assign a Gaussian prior distri- 
bution to 7, wc can preserve our least-squares approach to the fitting problem. We therefore fit 
by assuming that 7 = with an RMS uncertainty of Ubind = ih'md/V^- This value of cJbind has 
the same variance as the uniform distribution, and the entire ±7bind range (plus some marginally 
unbound orbits) is contained within the 2a contour, so we will obtain a realistic assessment of the 
uncertainties. 

Our procedure for handling the 7 degeneracy is thus: 

1. Fit the orbit with all six parameters free to vary. 

2. If the unconstrained fit yields S-y-y < cr^ind) then 7 is well constrained by the observations and 
we are done. Otherwise 7 is essentially decoupled from the observations, so we proceed using 
the binding constraint instead. 

3. Re-fit the observations but set 7 = and derive values and a covariance matrix for the 
remaining 5 parameters. Note that this best-fit orbit is very similar to the Vaisala orbit, as 
a line-of-sight velocity near zero implies an orbit near its perihelion or aphelion. 

4. Augment the covariance matrix by assigning E-y-y = cr^^^^, placing zeros in the remaining 
off-diagonal elements. 



The best-fit parameters a and covariance may then be used to predict positions and orbital 
elements as in the non-degenerate case. The 2a error ellipses will encompass > 95% of all bound 
orbits that arc consistent with the data, and not be restricted to circular or pcrihclic orbits. If 
we have some preconceptions about the maximal ellipticity of the orbits, this may be incorporated 
into the position predictions by placing a tighter a priori constraint on 7. 



2.5.2. Fitting Doubly Degenerate Orbits 

With only a single or few nights' data, or if there are only 2 observations, we are in the "slope 
regime." The above technique will fail due the additional degeneracy between a and 7 (and $ for 
observations out of the ecliptic plane). This will manifest itself in the covariance matrix obtained 
in Step 3 above, the 5-parameter fit to the orbit. The uncertainties T,aa and/or will be large. 

In this case we again introduce the binding constraint to limit the extent of the degeneracy. 
We express the transverse kinetic energy in terms of a parameter /;,: 

a'' + P^ = {l + h)GMe^\ (20) 

The binding constraint is then — 1 < < 1, with fb = corresponding to a nearly circular orbit. 
We crudely implement this constraint by minimizing a new quantity 

x' = x' + /6/3- (21) 

(with from Eq. 15), which pushes the solution toward a circular orbit, while yielding a covariance 
matrix that reflects both the observational uncertainties and the possibility of non-circular orbits. 
More sophisticated constrained optimizations are possible but there is little point to seeking precise 
constraints on such a poorly determined orbit. 



3. Implementation and Examples 
3.1. C Code 

All of the above algorithms have been implemented in a family of subroutines and drivers writ- 
ten in C. These inchidc routines for: creation and access of binary versions of the DE405 cphcmeris 
files (provided by David Hoffman); orbit integration; transformations among the various Cartesian 
and spherical coordinate systems; transformations between phase-space and orbital-element bases 
for orbits; reduction of positions to topocentric coordinates; and the fitting procedures described in 
the previous section. The code is portable and self-contained, and is available from the first author 
upon request. 

The main program is f it_radec, which takes as input a list of astrometric observations, and 
executes the orbit-fitting, including proper handling of degeneracies. The output is a file specifying 
the best-fit a description of the orbit, and its full covariance matrix E^. 
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A program abg_to_aei transforms the a parameters into traditional orbital elements. The 
covariance matrix is also transformed with a linear approximation. 

The CK information may be passed to the program predict, which generates predicted positions 
and their error ellipses for any specified observation date and site. 

The entire process is extremely fast on the typical present-day desktop computer. A fit to 
275 mock observations of Pluto over a 14-year arc takes about 1 second. The entire ensemble of 
observations of all 300 KBOs reported to the MPC are fit in a few seconds. 

3.2. Example and Verification 

We have verified the performance of the algorithms and software by conducting mock observa- 
tions of Pluto, Neptune, and several KBOs. Figure 2 shows the results of an ensemble of simulated 
observations of Pluto from 1995-1998. The true position of Pluto is obtained from the DE405 
ephemeris; mock observations are created at specified dates by adding Gaussian deviations to the 
true positions. An orbit is then fit to the mock observations. This fitted orbit is used to estimate 
the osculating barycentric elements a and e for Pluto as well as the uncertainty ellipse in the a—e 
plane. Alternatively we may use the orbit fit to estimate the (9x, Oy) sky position of Pluto at some 
future date and its projected error ellipse. We can then examine whether the estimated orbital 
elements or future position are consistent with the true ephemeris values. In all cases we find the 
mean value of the fit to the orbit to be as expected from the degrees of freedom in the fit. 

The first panel in Figure 2 shows the results of using 4 observations over an 11-day arc (JD 
2449927-2449938) to predict the location of Pluto fa 300 days after the first observation. Each 
point on the plot is the prediction from one of 1000 realizations of the mock 11-day arc. The 
star shows the true location of Pluto at the 300-day time. The ellipses show the la and 2a 
contours of the estimated error on the position given the observational errors. The estimated 
uncertainty ellipse is quite consistent with the scatter in the mock predictions. There appear to be 
two flaws in our prediction: first, the mean of the predicted positions does not coincide with the true 
position; second, the minor axis of the uncertainty ellipse appears to be significantly larger than the 
scatter. Both of these phenomena are expected, however, because the 11-day arc is the "acceleration 
regime" and is degenerate in 7. The fits therefore have assumed a value and uncertainty for 7. 
The displacement of the simulations from the true value is due to our assumption that 7 = for 
each fit, which is not the true value. But the uncertainty ellipse includes a contribution from our 
ignorance of 7, which makes the ellipse extend beyond the measured scatter, and in fact properly 
include the true position. 

In the second panel we simulate an observation sequence consisting of the 11-day arc plus 
and additional observation at 300 days. We plot the estimated a and e values from 1000 mock 
observations of the 300-day arc. In this case the uncertainty ellipses give a roughly correct estimate 
of the size of the errors, but the scatter does not match the error ellipses in detail. This demonstrates 
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Fig. 2. — Demonstration of position/orbital element estimation using the methods described in the text. 
In the left panels (a) and (c), mock observations of Pluto are used to predict a future position on the sky. 
Each small symbol is the prediction from one mock arc; the star is the correct position, while the two shaded 
ellipses are the 1- and 2-sigma uncertainty regions produced by our methods. The panel (b) shows the orbital 
parameters a and e estimated from 5 observations over 300 days, along with the true position and estimated 
uncertainties. Panel (d) shows the elements estimated after an additional observation at 1000 days. The 
scatter of estimates matches the uncertainties except in case (b). Details of this failure are explained in the 
text. 
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a limitation of our method, which is that the propagation of the covariance matrix Sq, to the orbital 

elements is incorrect if the map between them is significantly nonlinear. This is the case here, as 
Pluto is very near perihelion at this time. The 300-day arc determines the distance to Pluto quite 
well, and hence the scatter plot is bounded by a lower limit to the perihelion a(l — e). 

The third panel shows the results of using the 300-day arc to predict the position after 1000 
days, and the fourth panel shows the scatter in the a-e plane of elements determined from 6 
observations on a 1000-day arc. We can now usefully constrain all 6 parameters of the orbit, 
and the scatter about the true position and orbit are very accurately described by the estimated 
covariance matrices. 



4. Optimal Strategies 

An optimal schedule for recovery of KBOs would maximize some measure of the accuracy of 
the KBO orbit while minimizing the number of recovery observations required. A constraint on 
the observing schedule is that the positional uncertainty must not become larger than the field of 
view of the instrument, lest the recovery be missed and the object be lost. As an illustration of the 
utility of the methods described in this paper, we derive an optimal schedule of observations using 
a crude algorithm. 

We take our measure of the quality of the orbit fit to be the uncertainty cJa = V^aa in the 
semi-major axis of the orbit. This is not meant to be a universal metric for quality but rather a 
simple example. We place two constraints on the observing schedule: first, observations are not 
allowed when the solar elongation of the target is less than 90°. Second, we take the field of view of 
the recovery instrument to be ~ 10'; then the 1-sigma positional uncertainty must be < 2^5 if the 
full 2-sigma error ellipse is to be contained within the field of view for > 95% chance of recovery. 
Note that many CCD mosaics currently on 2-meter and 4-meter class telescopes have fields of view 
that substantially exceed 10'. 

A rigourous optimization for a set of N observations would consider all possible sets of N nights 
and choose the set which minimizes the final a a- We adopt a somewhat simpler "greedy" algorithm 
for optimization as follows: after observation M, we consider all nights for which observation M + 1 
would be possible, namely those satisfying the elongation ad positional uncertainty costraints. We 
next calculate the o"a which would result from an observation on each possible night; the M+1 
observation is then scheduled for the night which would result in minimal aa- 

Figure 3 illustrates the observation schedule derived for hypothetical discovery of Pluto in the 
year 1995. We assume that the object is discovered on a pair of nights roughly 45 days before 
opposition, then confirmed on another pair of nights 10 days later. We also presume that each 
observation has an error of cr = 0^'2. As seen in the Figure, the object must be observed once again 
in the discovery season as its uncertainty is too high to be recovered when it emerges from behind 
the Sun the following year. The object is now sufficiently well characterized to be found for nearly 
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3 years. When is the best time during this window to observe again? Interestingly the accuracy 

on a improves only weakly the longer we wait for our next observation. The "greedy" strategy 
suggests we wait 3 years, when the positional uncertainty approaches 2'.5. So after the initial 4 
discovery /confirmation images and only 2 recoveries, aa ~ 0.03 AU, and the orbit is known well 
enough to permit recovery for the next decade. 

A counterintuitive result is that the timing of further observations has little effect upon the 
accuracy of a — only increased number and/or accuracy of observations matters at this point. This 
is because the uncertainty in a derives primarily from the uncertainty in 7, i.e. the distance to 
the target. The leading observable consequence of 7 is the parallax ^yxEit), which does not grow 
with time beyond 1 year — hence the very weak improvement with longer arcs. Indeed we find that 
observing near the two quadratures of the second year yields a aa that is as low as spacing these 
two observations over several years. 

Another common observing scheme is to run discovery observations during two nights near 
opposition, then attempt confirmation one month later. For Pluto, the positional error is sufficiently 
small this second month to permit reliable confirmation. Henceforth the optimal strategy plays 
out just as the previous scenario: a further observation is required before quadrature to keep the 
object from being lost while behind the Sun, and the orbital quality depends primarily upon the 
number, not the arc length, of recoveries in later years. 

If the astrometric accuracy of the initial 11-day arc can be improved to C.'l, and the recovery 
tolerance raised to aa ~ 10', then no further observations are required in the first year. The object 
is recoverable at the beginning of the second season, two observations (near the quadratures) in 
this second year reduce aa to 0.0002 AU, and the position is good for a decade or more. 

These results are of course dependent upon the orbit of the object and upon the nature of one's 
metric for orbital quality. Our optimization method is, however, easily adapted to other metrics. 

5. Summary 

Orbit fits for KBOs and other distant objects are greatly simplified by the fact that the 
gravitational perturbation to the orbit is small compared both to the inertial motion and to the 
acceleration of the Earth. Gravity can be treated as a perturbation even for decades-long arcs, 
which simplifies the fitting process and the propagation of errors because there is no need to 
integrate the derivatives of the orbit. The orbit itself can be integrated quickly because large time 
steps can be taken. Distant objects move slowly across the sky, so we can also use the convenience 
of a tangent-plane projection of the celestial sphere. The dominance of reflex motion over proper 
motion means that KBO orbit solutions are unambiguous. 

Exploiting these properties, we have created algorithms and software which quickly and ac- 
curately calculates orbital elements and ephemerides and their associated uncertainties for targets 
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Fig. 3. — The observation optimization routine described in the text is apphed here to Pluto. We assume 
an initial set of four observations over a 10 days, about 6 weeks before opposition. We then ask when further 
observations are best scheduled. The upward-sloping lines depict the uncertainty in position of the object 
given the 10-day arc; another observation must be done while the uncertainty is below 150" (the dotted line) 
in order to assure successful recovery in a 10' FOV. The curved is dashed when the target is too close to 
the Sun to observe, solid otherwise. The downward-sloping lines show the uncertainty in semi-major axis 
a that would be obtained from a further observation at the given time. Additional observations (stars) are 
scheduled at the time during the recovery window when the measurement will most reduce the uncertainty on 
a. Each new observation then starts a new pair of curves for a and position uncertainty. A new observation 
is needed at 4 months (before the target disappears behind the Sun), after which the uncertainty on a is 
only weakly dependent upon the timing of observations. Another observation is needed within 3 years else 
the target becomes lost. 
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> 10 AU from the Sun. In cases where the observed arc is short enough to leave the orbit de- 
generate, we can still calculate sensible positional errors because wc have chosen an orbital basis 
in which the degeneracy is confined to one parameter, namely the line-of-sight velocity. For very 
short arcs in the "slope regime," 2 additional parameters become degenerate. We are still able to 
place error ellipses on ephemerides which are conservative in the sense that they will include any 
bound orbit. Most short-arc objects will be recovered closer to the MPC predicted ephemeris than 
our error ellipses would suggest, because the philosophy of the MPC is to use the stronger a priori 
consideration that the orbits resemble those of the majority of the known populations. 

Our routines should prove valuable for planning of recovery observations, as we have illustrated 
with some simple examples. Estimates of the uncertainty ellipsoids in orbital-element space will 
also be needed for study of the long-term dynamics of the known KBOs. 

In the near future we expect that KBO positions will be measured to accuracies of a few milli- 
arcseconds with orbiting observatories — though such observations will be even more expensive than 
present-day KBO observations. It will then be possible to place strong constraints on KBO orbits 
even from very short {e.g. 24-hour) arcs. For example, the motion of an observatory in low-Earth 
orbit imparts a reflex motion on the KBO that is easily detected at this level, giving a measure 
of the distance to the target in less than an hour. The methods described herein will provide the 
means to plan and to analyze such observations. 

We thank David Hoffman for making available his C code to read the DE405 ephemerides; Brian 
Marsden for discussions on fitting of KBO orbits; Julia Plummer for assistance with Figure 1; and 
Remi Malhotra and Lynne Allen for serving as "guinea pigs" for tests of the software. This work 
is supported by NASA Planetary Astronomy grant #NAG5-7860, and grant #AST-9624592 from 
the National Science Foundation. 
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A. Coordinate Transformations 

A.l. Angular Coordinates 

We make use of three sets of coordinates on the celestial sphere: equatorial coordinates (a, 5) 
in the ICRS reference frame, J2000; ecliptic latitude and longitude {l,b), also taken in the epoch 
2000 frame; and our projected coordinates {9x,0y) measured in a tangent plane about some ref- 
erence direction (q!0)<^o) or (£0,^*0)- The transformation from equatorial to ecliptic coordinates is 
implemented as 

sin b = cos e sin 6 — sin e cos 5 sin a (Al) 

cos e cos 6 sin a + sin e sin 6 , , „, 

tan^ = , (A2) 

cos cos a 

where e = 23? 43928 is the obliquity of the ecliptic at J2000 (Standish et al. 1995). The inverse 
transformation simply reverses the sign of e. 

The map from ecliptic to tangent-plane coordinates about (^O; ^0) is 

^ ^ cos b sin ji - ip) 

^ sin 60 sin 6 — cos 60 cos 6 cos (£ — ^o) 



cos bo cos b — sin 60 sin b cos {i — io) 
sin bo sin b + cos bo cos b cos {£ — io)' 



(A4) 



. sin 60 + 6*2/ cos 60 /.^s 
smO = , =- (A5) 



The inverse transformation is 



1 + 02 + ^2 

■re e \ ^x/ cos 60 

1 + ^2 + ^2 

Partial derivative matrices are straightforward for all transformations, and approximations for 
Ox-,dy^ 1 are particularly simple. 



A. 2. Spatial Coordinates 

There are also three relevant spatial coordinate systems. The equatorial Xcq has origin at 
the Solar System barycenter, with x positive to the Vernal Equinox and z positive to the North 
equatorial pole of J2000. The ecliptic system also has origin at the barycenter, but z is positive 
to the North ecliptic pole of J2000. Our "telescope-centric" coordinate system has origin at the 
location of the observer at t = 0, with z positive along the line of sight to the target at t = 0, the 
y axis in the plane of z and the North ecliptic pole — see Figure 1. 

The transformation from equatorial to ecliptic coordinates is 
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The inverse transformation inverts the sign of e. 

The transformation from echptic to telescope-centric coordinates is 
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The ecUptic-coordinate location of the observer at time t = is given by the vector {a;o,yO)^o}- 
The inverse transformation is 
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Since all the transformations are linear, the partial derivatives are trivial. 



(Bl) 



B. Orbit Basis Transformations 

The orbit-fitting algorithms produce a description of the orbit in terms of ct = {a, /3, 7, d, /3, 7}. 
To obtain orbital elements, we first transform to Cartesian phase space elements in the telescope- 
centric coordinate system: 

XT = 01/^ ut = (9/7 zt = 1/7 

XT = a/^ jjT = $/j Zt = 7/7 
The telescope-centric phase space coordinates can then be transformed to barycentric ecliptic Carte- 
sian phase-space coordiates using Eq. (A9). The transformation matrix for the velocities is the same 
as that for the position, but there is no telescope-to-barycenter ofi^set for the velocities. 

With the barycentric position and velocity x and v in hand, the osculating orbital elements 
are determined with the standard formulae below. Using the gravitational constant fj, = GM, we 
have 



Semimajor axis: 



Eccentricity: 
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With the further definitions 



h = XXV (B5) 
n = z X h (B6) 



the inchnation, ascending node, and argument of periheUon are given by 



cosi = ^ (B7) 

71 

cosn = — (B8) 
n 

n • e 

cos a; = . (B9) 

ne 



The eccentric and mean anomaUes are 



t&nE = -4^ — (BIO) 
M = E-esinE, (Bll) 



with 



p — X 



e 



Finally, the time of periapse passage is 



(B12) 



y = — / (B13) 

b = a\J\ - e2 (B14) 
p = — (B15) 



Tp = tQ- M^^fi/Ji, (B16) 

where to is the time at which x and v are determined. 

The partial derivatives necessary for propagation of errors are tedious but calculable. Our 
phase-space elements are centered on the solar system barycenter, therefore the osculating elements 
are also barycentric. We take M as total mass of the solar system since we are using the barycenter. 
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